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SOME  EXPERIMENTS  IN  ALGEBRAIC 
MANIPULATION  BY  COMPUTER 

ABSTRACT 

A  set  of  subroutines  to  allow  algebraic  manipulations 
on  the  IBM  7094  computer  has  been  written  using  a  List 
Processor,  SLIP. 

A  series  of  four  problems  of  Increasing  difficulty 
were  solved  using  these  routines.   These  problems  were: 

1)  To  expand  and  simplify  a  polynomial  in 
three  variables . 

2)  To  substitute  one  truncated  power  series 
into  another. 

3)  To  obtain  the  series  solutions  of  two  ordinary 
differential  equations  by  equating  coefficients . 

4)  To  obtain  a  partial  expansion  of  the  Fourier  series 
solution  to  Laplace's  equation  on  a  rectangular 
region  with  one  irregular  boundary. 
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NYO-li^  80-11 
SOME  EXPERIMENTS  IN  ALGEBRAIC  MANIPULATION  BY  COMPUTER 

by 

Arnold  Lapldus,   Max  Goldstein 
Introduction. 


Symbolic  algebraic  manipulation  is  an  area  of 
non-numerical  application  of  digital  computers  which  has 
recently  attracted  considerable  interest.   In  fact, 
several  groups  have  already  written  programs  and  systems 
of  this  type.   Two  highly  developed  systems  are  Bell's 
ALPAK  [1]  and  IBM's  FORMAC  [2].   ALPAK  does  non-numerical 
algebra  on  polynomials  with  an  arbitrary  number  of 
variables.   FORMAC  is  an  unreleased  experimental  formula 
manipulation  compiler.   Since  neither  of  these  systems 
was  directly  applicable  to  our  needs,  and  to  gain 
experience,  we  decided  to  develop  our  own  routines  for 
algebraic  manipulation.   We  call  the  collection  of  these 
routines  for  the  IBM  709^,  SYMB0LANG . 

Limited  work  had  already  been  done  here  in  the 
non-numerical  area.   Successful  fluid  dynamics  computa- 
tions using  analytic  continuation  [^,^,5]    led  ultimately 
to  a  package  of  routines  for  manipulating  power  series 
in  two  variables  [6].   A  symbolic  differentiation  subroutine 
had  also  been  written  [7] •   Neither  of  these  was  intended 
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to  be  a  general  purpose  manipulator;  the  power  series 
package  did  not  allow  the  use  of  explicit  variables  and 
the  differentiator  did  not  allow  further  algebraic  or 
numerical  manipulation. 

The  strategy  we  adopted  for  this  project  was  to 
attempt  to  solve  a  sequence  of  algebraic  problems  of 
increased  difficulty.   The  necessary  subroutines  would 
then  form  the  basis  for  further  work  in  developing  a 
capacity  for  algebraic  manipulation  applicable  to  an 
increasing  variety  of  problems.   Fiarthermore,  this 
system  was  to  be  experimental.   We  realized  that  as 
soon  as  a  few  problems  were  solved  we  would  probably 
want  to  use  the  information  obtained  from  a  first  try 
to  write  a  more  efficient  and  useful  set  of  routines. 

At  the  outset,  the  points  which  needed  considera- 
tion were  those  which  dealt  with  the  representation  of 
algebraic  expressions;  in  particular,  notation  and 
memory  organization.   These  points  are  discussed  in 
Sections  2  and  5» 

Section  4  consists  of  descriptions  of  routines 
which  were  developed  for  algebraic  manipulation.   With 
these  routines,  four  problems  were  attempted  and  solved. 
These  problems,  in  order  of  difficulty,  were: 

1)  To  expand  and  simplify  a  polynomial  in  three  variables 

2)  To  substitute  one  truncated  power  series  into  another. 

3)  To  obtain  the  series  solutions  of  two  ordinary 
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differential  equatlonsby  equating  coefficients. 
4)  To  obtain  a  partial  expansion  of  the  Fourier  series 
solution  to  Laplace's  equation  on  a  rectangular 
region  with  one  irregular  boundary. 

Their  solutions  are  described  in  Section  5. 

2.   Algebraic  Expressions  -   Notation,  List  Processors. 

A.  Notation 

In  devising  a  notation  for  algebraic  expressions, 
both  the  internal  and  external  forms  needed  consideration. 
Obviously,  an  external  notation  which  could  be  easily 
written  and  read  by  the  programmer  was  needed,  and  because 
most  computer  references  to  an  expression  would  be  in  the 
processes  of  m^ultiplication  and  collection  of  terms,  an 
internal  representation  which  facilitated  the  execution 
of  these  basic  processes  was  desirable.   In  order  to 
simplify  debugging  in  this  first  attempt,  it  seemed  wise 
to  use  a  representation  which  was  essentially  the  same 
internally  and  externally.   This  would  have  two  other 
advantages  -■=  complicated  conversion  routines  would  not 
have  to  be  written  and  there  would  be  sufficient  flexi^ 
bllity  available  in  the  notation  so  that  it  could  be 
expanded  to  Include  notation  for  a  larger  class  of 
variables.   For  example,  in  the  first  problems  discussed 
in  Section  5,  only  algebraic  variables  such  as  x,  y,  ... 


had  to  be  processed,  but  in  the  last  problem,  cosine 

combinations  appeared  in  the  expressions. 

The  notation  which  was  actually  used  in  SYMB0LANG 

is  the  following: 

Variable:  A  variable  such  as  x,y, ...  is  denoted  by  a  full 
word  of  Hollerith  data  which  starts  with 
two  blanks. 

Variable  Raised  to  a  Power:   If  a  variable  is  raised  to 
a  power  k  7^  1,  then  the  variable  is  followed 
by  a  full  Hollerith  word  bbbb-^^,  and  then 
another  full  word  which  contains  k  in  floating 
point.   If  the  variable  is  raised  to  the  power  1, 
then  it  is  not  followed  by  the  Hollerith  word 
bbbb**,  nor  is  it  followed  by  1;  i.e.,  the 
exponent  is  understood  to  be  1  (b  represents 
blank) . 

Coefficient:   A  coefficient  is  a  normalized,  non-zero, 
floating  point  number. 

Term:     A  term  is  a  coefficient  followed  by  variables 
raised  to  powers . 

Expression:   An  expression  is  a  sequence  of  terms  on  a 
single  list. 

Example;   4yx  -  6y  x  "^  would  appear  on  a  list  as 
follows: 
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Mixed  Format 


Octal 


4.0 

203400000000 

bbbbbY 

606060606070 

bbbbbX 

606060606067 

bbbb** 

606060605454 

2. 

202400000000 

-6. 

605600000000 

bbbbbY 

606060606070 

bbbb** 

606060605454 

2. 

202400000000 

bbbbbX 

606060606067 

bbbb** 

606060605454 

-3. 

602600000000 

Notice  that  exponents  as  well  as  coefficients  are 
floating  point  numbers,  so  there  may  appear  to  be  some 
ambiguity  if  the  list  is  glanced  at  superficially.   How- 
ever, if  the  list  is  read  in  sequence  the  apparent 
ambiguity  of  the  floating  point  numbers  is  resolved. 
Note  also  that  exponents  may  be  negative  as  well  as 
positive . 

In  this  notation,  no  provision  is  made  for  factored 
expressions.   The  reason  for  this  is  that  subroutines 
which  could  manipulate  such  expressions,  as  well  as  those 
which  would  be  in  expanded  form,  would  be  complex  and 
time  consum.ing.   So  far,  a  factorization  capability  has 


not  been  necessary  because  the  end  results  of  the  computa- 
tions have  been  tractable  expressions  in  their  expanded  form. 

We  now  see  that  a  finite  polynomial  expression  can  be 
represented  conveniently  by  a  string  of  computer  words. 
A  convenient  computer  construct  for  dealing  with  such 
strings  is  available  in  what  will  here  be  called  LISTS, 
i.e.  the  lists  used  in  list  processors. 

B.   List  Processor 

The  decision  to  use  a  list  processor  in  this  experi- 
mental program  was  based  exclusively  on  the  fact  that  list 
processors  allow  an  extremely  flexible  dynamic  memory 
organization. 

In  order  to  develop  the  concept  of  LIST'S  as  used 
in  list  processors,  they  will  be  compared  to  Fortran  arrays. 
While  some  of  the  advantages  of  LIST  operations  can  be 
obtained  from  Fortran  arrays  [12],  list  processors  allow 
these  operations  to  be  implemented  very  simply;  e.g.,  as 
a  CALL  or  a  r([ACR0  or  as  a  subexpression  in  the  right  side 
of  a  Fortran  statement. 

The  list's  which  are  the  objects  of  list  processors 
are  similar  to  Fortran  arrays  in  that  each  element  of  a 
LIST  usually  contains  a  single  word  of  data. 

It  is,  of  course,  the  differences  between  LIST'S  and 
arrays  which  are  im.portant .   Some  of  them  follow: 
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Property  1.   There's  always  room  for  one  more  element  on 
a  LIST,  provided  there's  room  for  an  element 
on  any  LIST  In  the  memory.   No  LIST  is  ever 
too  long  unless  the  whole  memory  set  aside 
for  LIST'S  is  full  with  current  information. 

Property  2.   A  LIST  can  be  created  at  any  time  in  the 
running  of  the  program.   Distinguish  this 
property  from  an  array  which  must  be  set  up 
prior  to  compiling. 

Property  3-   When  a  LIST  is  no  longer  needed,  and  destroyed, 
no  trace  of  it  is  left  (except  possibly  for 
debugging  purposes Ji,  and  it  no  longer  takes 
space . 

Property  4.   When  an  element  on  a  LIST  isn't  needed  any 
more  and  is  removed  from  a  LIST,  it,  too, 
no  longer  takes  space. 

Property  5.   Any  LIST  may  be  converted  into  a  pushdown 
list*  or  stack. 

For  an  algebraic  manipulator  in  which  the  number  of 
terms  appearing  in  intermediate  expressions  is  large  and 
unknown.  Property  1  is  indispensable.   Properties  3,  4,  5 
allow  subroutines  to  create  intermediate  algebraic 
expressions  of  arbitrary  size  without  using  any  part  of 


* 
Pushdown  List:  A  list  of  elements  such  that  a  new 

element  is  put  on  top  and  when  an  element  is  retrieved 

from  the  list,  it  is  the  top  element;  i.e.,  LAST  IN/ 

FIRST  OUT. 
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the  memory  permanently,  and  without  setting  aside  a  working 
storage  region  which  would  remain  unused. 

These  properties  are  available  in  the  list  processor 
called  SLIP  [8].   The  reasons  for  choosing  SLIP  instead 
of  another  list  processor  were: 

1.  Its  list  processing  facilities  are  as  varied  and 
appropriate  for  our  needs  as  any  other  list  processor. 

2.  It  is  superbly  documented. 

3.  It  is  written  as  a  collection  of  Fortran  subroutines. 

The  Fortran  organization  of  SLIP  also  allows  subroutine 
additions  to  SLIP  without  interfering  with  previously  written 
routines.   This  capability  was  used  in  designing  an  output 
routine,  and  in  changing  some  of  the  SLIP  routines . 

The  fact  that  SLIP  was  written  in  Fortran  allowed  the 
algebraic  manipulator  to  be  written  in  Fortran,  and  it 
therefore  has  the  same  advantages  that  other  types  of 
Fortran  programs  have.   Unfortunately,  it  has  the  same 
major  disadvantage  but  to  a  higher  degree  =-  it  takes  much 
more  running  time  than  would  a  FAP  designed  program.   The 
reason  for  this  is  that  each  operation  on  a  LIST  requires 
many  references  to  whole  Fortran  subroutines  with  their 
index  storing  and  address  setting  paraphernalia,  and 
often,  the  SLIP  system  goes  through  as  many  as  ten 
subroutine  levels  to  perform  an  operation  on  a  LIST. 

We  traded  running  time  for  convenience,  flexibility 
and  clarity  in  programming;  we  used  SLIP  for  a  list 
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processor  and  coded  the  manipulator  as  Fortran  subroutines 

3.   SLIP  Lists  In  SYMB0LANG. 

The  objects  upon  which  the  algebraic  manipulator 
operates  are  algebraic  expressions .   As  explained 
previously,  these  appear  on  LIST'S  of  the  SLIP  list 
processor.   In  this  section  portions  of  the  SLIP  system 
which  are  needed  or  useful  In  coding  a  problem  will  be 
explained.   But  first  a  discussion  of  LIST  names  Is 
necessary. 

While  the  list  processor  doesn't  actually  store 
successive  elements  of  a  list  In  sequential  locations 
In  the  computer,  a  list  will  be  assumed  to  have  the 
form  of  an  ordered  sequence  of  computer  words. 

A.   LIST  Names 

References  to  a  LIST  must  be  made  via  a  LIST  NAME. 
The  NAME  of  a  LIST  Is  first  created  at  the  same  time 
that  the  LIST  Is  created.   A  NAME  Is  a  pseudo-Integer; 
It  must  follow  the  same  rules  for  naming  that  a  Fortran 
Integer  does,  but  It  cannot  be  used  In  arithmetic  state- 
ments except  that  any  Fortran  Integer  can  be  set  equal  to 
a  previously  defined  NAME  and  It  then  becomes  a  NAME  also. 

Example  1.   The  FORTRAN  statement: 

LIST  1  =  LIST  (9) 
causes  the  subroutine  LIST  (9)  to  assign  the  Fortran 
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variable  LIST  1  the  value  of  a  LIST  NAME.   LIST  (9)  also 
sets  up  a  LIST. 

Example  2.   The  FORTRAN  statement: 

L3  =  LIST  1 
causes  the  Fortran  variable  L3  to  become  a  NAME  referring 
to  the  same  LIST  as  the  already  defined  LIST  1.   NOTE:   If 
two  NAMES  refer  to  the  same  LIST  and  one  of  them  is  changed 
to  refer  to  another  LIST,  then  the  second  of  them  will 
continue  to  refer  to  the  first  LIST. 

B,   SLIP  Subroutines  for  Creating,  Filling,  Emptying, 
Destroying  and  Testing  of  LIST'S. 

The  following  are  descriptions  of  som.e  SLIP  subroutines 

which  are  directly  applicable  to  the  LIST's  that  are  used  in 

algebraic  m.anipulation: 

i.    LIST  (9) -The  statement  LI  =  LIST  (9)  causes  the 

function  LIST  (9)  to  create  an  empty  LIST  and  leave 
its  NAME  in  the  Fortran  variable  LI.  This  is  the  only 
function  of  this  collection  which  has  a  dummy  variable. 

li.   LSSCPY  (LI)  -  The  input  parameter  LI  must  be  a  NAME. 
The  function  creates  a  LIST,  the  NAME  of  which  is 
the  value  of  the  function.   The  newly  created  LIST 
is  a  copy  of  the  input  LIST  LI. 

iii.  MANY(L1,A,B,C  ...  ,N,M,  ..<.)- The  input  parameter  LI  must 
be  a  NAME.   The  variables  A,B, . . .  are  arbitrary  in 
number  and  are  Fortran  variables.   MANY  puts  the 
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variables  A,B, . . ,  onto  the  bottom  of  the  LIST  LI  In 
the  order  that  they  are  listed. 

Iv.   MTLIST  (NAME)-  The  LIST  the  NAME  of  which  appears  as 
an  Input  parameter  to  this  function  Is  emptied  --=  Its 
memory  space  Is  made  available  for  other  LIST'S.   The 
value  of  the  function  Is  the  NAME  of  the  LIST. 

V.    IRALST  (NAME)  -  The  LIST  whose  NAME  appears  as  an 
Input  parameter  Is  erased.   All  the  memory  which 
It  was  using  Is  made  available  to  other  LIST'S. 
The  value  of  the  function  Is  zero. 

vl.   LSTEQL  (L1,L2)  -  The  two  Input  parameters  are  NAME'S 
of  LIST'S.   If  the  two  LIST'S  contain  the  same 
elements  In  the  same  order,  the  value  of  the  function 
Is  zero. 

vli.  LISTMT  (LI)  -The  Input  parameter  LI  Is  a  NAME.   If 
It  Is  the  NAME  of  a  LIST  which  Is  empty  then  the 
function  value  Is  zero;  otherwise.  It  Is  non-zero. 

C.   Output  of  Lists 

The  contents  of  a  LIST  Is  obtained  by  a  statement 
of  the  form: 

CALL   OTJTLWR  (  6Hxxxxxx,L1) 
where  any  six  Hollerith  characters  can  be  substituted 
for  xxxxxx  and  LI  Is  the  NAME  of  the  LIST  to  be  printed. 
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The  statement  causes  the  computer  to  print  output 
similar  to: 

075211075211 
This  LIST  is  XXXXXX 

-1.00000000     6oi4oooooooo     060107000000 
1.00000000     201400000000     066541000000 

A  606060606021      056171000000 

In  the  above,  the  first  number  printed  is  the  Fortran 
variable  corresponding  to  the  NAME  of  the  list  LI.   The 
word  XXXXXX  is  the  word  of  Hollerith  data  supplied  to  the 
output  routine.   All  the  data  in  the  LIST  is  then  printed 
as  follows:   Position  nine  of  the  data  word  is  examined j 
if  it  is  a  zero  the  data  is  assumed  to  be  Hollerith  and 
printed  as  6  Hollerith  characters,  if  the  bit  is  a  1, 
then  the  data  is  printed  as  a  floating  point  number.   In 
both  cases  the  octal  representation  of  the  word  and  the 
location  of  the  SLIP  cell  containing  the  data  in  709^ 
memory  are  printed  on  the  same  line. 

4.   Manipulative  Routines » 

The  manipulative  routines  which  have  been  written  as 
part  of  SYMB0LANG  fall  into  two  major  classes: 

(1)  those  which  do  manipulations  usually  recognized  as 
algebraic  —  e.g.,  multiplication,  substitution,  and 

(2)  those  programmed  to  do  algebraic  operations  which 
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are  so  simple  In  pencil  and  paper  calculations  that  they 
usually  wouldn't  he  called  algebraic  •=-  e.g.,  copying  a 
term,  finding  the  exponent  of  a  variable  In  a  term, 
separating  an  expression  Into  terms  which  do  contain  a 
specific  variable  and  those  which  do  not,  etc.   In  this 
section,  the  calling  sequences  for  Fortran  subroutines 
and  functions  which  do  operations  of  the  algebraic  type 
are  explained  first,  and  then  those  of  the  second  type 
are  explained. 

A.   Manipulative  Routines  of  the  Algebraic  Type. 

1.   TRMMPY  (LSTA,LSTB,LSTC)  -  Assumes  LSTA,  LSTB  are 
both  names  of  LIST'S  which  contain  one  term 
each.   The  subroutine  multiplies  the  terms 
and  puts  the  result  onto  the  LIST  which  has 
the  NAME  LSTC .   The  LIST'S  LSTA,  LSTB  are  not 
changed . 

11.   SUMPY  (LSTA,  LSTB,  LSTC)  -  LSTA,  LSTB  are  both 

names  of  LIST'S  which  contain  expressions. 
This  subroutine  multiplies  the  expressions 
on  LSTA  and  LSTB  and  adds  the  results  onto 
LSTC.   The  contents  of  the  LIST'S  LSTA, 
LSTB  are  not  changed , 

ill.  DVTRM  (LSTl,  VRBL)  ~  LSTl  is  a  NAME  of  a  LIST  which 
contains  exactly  one  term.   VRBL  is  to  be  a 
floating  point  number  of  a  Hollerith  variable, 
This  subroutine  divides  the  term  on  LSTl  by 

VRBL  and  puts  the  result  back  on  LSTl. 
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iv.   DVSUM  (LST1,VRBL,LST2)  -  LSTl,  LST2  are  NAME'S  of  LIST'S. 
VRBL  is  either  a  floating  point  variable  or  a 
Hollerith  variable.   This  subroutine  divides 
the  expression  on  LSTl  by  VRBL  and  adds  the 
results  to  LST2.   The  expression  on  LSTl  is 
not  altered. 

V.    SBST  ( LSTl, CHAR, LIST2, LIST3 )~  LISTl,  LIST2,  LIST5 

are  LIST  NAME'S.   CHAR  is  a  Hollerith  variable. 
The  subroutine  substitutes  the  expression  on 
LISTl  for  the  variable  CHAR  wherever  it  appears 

-  -   -'    on  LIST2  and  adds  the  result  to  LIST3.   The 
expression  on  LISTl  is  not  changed. 
Cautionary  note:   the  expression  on  LISTl  should 
not  contain  the  character  CHAR. 

vi.   DERIV  (LISTl, CHAR, LIST2)  =  This  subroutine  differenti- 
ates the  polynomial  on  LISTl  with  respect  to  the 
Hollerith  variable  CHAR  and  adds  the  result  to 
LIST2,   The  polynomial  on  LISTl  is  left  unchanged 

vii.  GETC0E  ( CHAR, LEXP, LISTl, LIST2)  -   CHAR  is  a 

Hollerith  variable.   LEXP  is  a  non-negative 
Fortran  integer.   LISTl  is  the  NAI'-E  of  a  LIST. 
LIST2  is  the  NAME  of  another  LIST.   The 
subroutine  obtains  the  coefficient  of  CHAR  ** 
LEXP  on  LISTl  and  adds  it  to  LIST2.   LISTl  is 
left  unchanged . 
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viil.SMPL(LSTC)  -  LSTC  Is  a  NAME.   The  subroutine  SMPL 

collects  terms  on  LSTC  and  puts  the  simplified 
expression  back  onto  the  list. 

Ix.   TRUNC(L1,VAR,PWR)  -   Ll  Is  a  LIST  NAME,  VAR  a 

Hollerith  variable,  PWR  a  Fortran  floating 
point  number.   The  subroutine  erases  all 
terms  on  Ll  that  contain  the  variable  VAR 
to  a  power  which  is  greater  than  or  equal  to 
PWR. 

X.    TOCAL  (LIST,T,C0EP)  -  LIST  is  a  LIST  NAME.   T  and 

COEP  are  Hollerith  variables.   This  subroutine 
searches  LIST  for  the  smallest  exponent  M 
such  that  C0EP  and  T-**-*M  are  contained  in  the 
same  term.   It  then  deletes  each  term  on  LIST 
which  contains  T  to  a  power  larger  than  M. 

xl.   SOLVE  (LIST1,CHAR,LIST2)  -  LISTl,  LIST2  are  LIST 
name's;  char  is  a  Hollerith  variable.   For 
this  subroutine,  it  is  assumed  that 

1.  the  expression  on  LISTl  is  equal  to  zero. 

2.  CHAR  appears  linearly  and  in  only  one 
term  of  the  simplified  expression  on  LISTl, 

This  subroutine  solves  for  CHAR  and  puts  the 
result  onto  LIST2,   LISTl  is  not  altered. 


B.  Manipulative  Subroutines  of  non-algebraic  type. 

i.    BREAK  (LST1,CHAR,LST2,LST3)  "  LSTl,  LST2,  LST3  are 

all  LIST  NAME'S.   CHAR  is  a  Hollerith  variable, 
This  subroutine  separates  the  terms  of  LSTl 
into  those  which  contain  the  variable  CHAR 
and  those  which  do  not.   Those  terms  which 
have  CHAR  in  it  are  put  on  LST2,  those  which 
don*t  are  put  on  LST3.   LSTl  is  left  as  it 
was  originally. 

ii.   CPYTRM  (LSTPRM,LSTT0)  -  This  function  gets  the 
value  0  if  LSTPRM  was  originally  empty  and 
the  value  1  otherwise.   If  LSTERM  isn't 
empty,  then  a  term  is  copied  from  the  top 
of  LSTMIM  to  the  bottom  of  LSTT0. 

iii.  P0l^ffiR  (Ll.VAR)  -  LI  is  a  LIST  NAME,  VAR  is  a 

Hollerith  variable.   If  the  first  term  on  LI 
which  contains  the  variable  VAR  contains  it 
to  the  power  P,  then  the  function  value 
assigned  to  P0WER  is  P .   If  VAR  doesn't 
appear  on  the  LIST  LI,  then  the  function 
value  is  zero.   LI  is  left  as  it  was.   This 
function  is  usually  used  in  the  case  that  LI 
has  only  one  term  on  it. 

iv.   SRIRM  (LSTl)  -  LSTl  is  a  LIST  NAME.    LSTl  contains 
exactly  one  term.   This  subroutine  arranges 
the  variables  of  a  term  in  alphabetic  order. 
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C.   A  Very  Simple  Example 

Complete  coding  to  obtain  (A+B+C)   in  expanded 
form  follows: 

C0MM0N    AVSL,W,X 

DIMENSI0N  W(IOO),  X(IOOO) 

CALL  INITAS  (X.IOOO) 

1  LI  =  LIST(9) 

2  L2  =  LIST(9) 
5     L3  =  LIST(9) 

4  CALL  MANY  (LI,  1. ,  6HblDbbbD,  6Hbbbb**,5.) 

5  CALL  MANY  'L2, 1 . , 6HbbbbbA, 1 . , 6HbbbbbB, 1 . , 6HbbbbbC) 

6  CALL  SBST  (L2,6HbbbbbD,Ll,L3) 

7  CALL  0UTLWR  (6HANSWER,L3) 
CALL  EXIT 

END 
The  common  and  dimension  statements  are  required  by 
the  SLIP  system.   Note  that  the  Fortran  variables  W  and  X 
are  not  available  to  the  programmer.   The  CALL  INITAS 
statement  is  also  required  by  the  SLIP  system  for  initiali- 
zation.  Statements  niimbered  1,  2,  3  set  up  three  empty 
LIST'S  which  have  the  NAME'S  LI,  L2,  L3-   Statement  4 
puts  the  expression  D   onto  LI.   Statement  5  puts 
the  expression  A+B+C  onto  L2.   Statement  6  substitutes 
A+B+C  for  D  and  puts  the  result  onto  the  LIST  that  has 
the  NAME  L3.   Statement  7  prints  out  the  result. 
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5.   Problems  Solved  Using  the  Algebraic  Manipulator. 

This  section  consists  of  a  description  of  problems 
solved  using  the  subroutines  which  were  described  in  the 
previous  section.   They  are  arranged  in  order  of  increas- 
ing difficulty.   The  second  problem  is  harder  than  the 
first  because  it  requires  truncation.   The  third  problem 
uses  an  implicit  type  of  truncation  and  also  solves  a 
sequence  of  simple  linear  equations  as  part  of  the 
problem.   In  the  fourth  problem  additional  routines  had 
to  be  developed  which  could  deal  with  elementary  combina- 
tions of  the  cosine  functions. 

Ao   Polynomial  Manipulation. 

This  problem  occurred  as  part  of  the  stability 
analysis  of  a  difference  scheme  [10]. 

Let  A  =  4v^(A^+l-2A+2Av^-v^)[A^(l~v)^+2A(l-v)(l+v^)(v-l-v^) 
+  (l+v^)^(v-l-v^)^] 

Let  B  =  v|j(v^+l)^[2A-2+v(5-v)]^+4v2(2-v)^v^(v^+l)^(A+v^-l)^ 
+  (1-v)^v^(2A^+v|J+2-4a+4Av^-3v^)^-4v^(v^+1)^(2-v)- V 
(2A-2+3v-v^)(A+v^-l)+2v^(v^+l)^(l-v)v(2A-2+3v-v^) • 
(2A^+v,V2-4a+4Av^=5v^)=4v^(2-v)  (l-v)v^(v^+l)  • 
(A+v^-1)(2a2+vJ^+2-4a+4Av^-3v^)  . 

¥e  wanted  to  form  the  expression  B-A  and  check  the 
coefficients  of  7\,    I.e., 
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(1)  The  coefficient  of 

(2)  The  coefficient  of 
1?   =  8v^v(l+v^)^(l"v) 

(3)  The  coefficient  of 

(4)  The  coefficient  of 

A  =  4v^fl+v^}^(-2v^-4v^»2vU"12v^-8v^v^-2v^+8v^+2v^v^ 

+  6v+8v,  v+4v,  v-i-2v?v) 
h    h    h 

(5)  The  coefficient  of 

A°  =  4v^^(l+v^)^(vh+2v^+v|j+v^-4v^-2v^v^+5v^+6v^v^+2v^v^ 

2      ^ 
-  2v  ~  4v,  V  -  4v,  V  -  2vr,v) 
h      h      h 

Notice  that  the  expressions  (1=  v),  (1+v,  ),  (v^l-v,  ), 
(A+v^-1),  (2A™2+3v-v^),  (2-v)  and  (2A^+v|J+2-4A+4Av^"3v^) 
all  appear  more  than  once.   These  were  replaced  by  simple 
variables  in  expressions  A  and  B  so  that  substitutions 
could  be  made  using  the  SBST  routine.   Expression  B-A 
was  formed,  expanded  and  simplified  and  then  the  coefficients 
of  the  powers  of  A  were  printed  out.   The  coefficients  (l)-(5) 
were  expanded  by  hand  and  compared  with  the  computed  coeffici- 
ents.  They  were  correct. 
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B.   Truncated  Power  Series  Substitution. 

Let   e^^^    =   1   +  z   +  z^/21    +   ...    +  z^/\cl 

Let   sin  y  =   y  -  y^/3l    +  y^/5J    -    •••    +  yVkl 

The  problem  was  to  find  e^       by  substituting 

(z) 

the  polynomial  for  sin  y  into  the  polynomial  for  e^   , 

The  result  was 

g(sin  y)  ^  1  +  y  ^  ^5y2  __  ^^25y^  _  ^q^^^^^^^i^^^ 


-   .oo4l666663y  +  .oiiiiiiioy^  . 


Because  the  computation  was  performed  using  single 
precision  floating  point  arithmetic  on  the  ^Q^k ,    the 
eighth  significant  figure  is  not  accurate. 

C.   Bessel  Functions  for  Large  Arguments. 

The  following  mathematical  analysis  is  a  summary  of 

portions  of  [11] . 

Let   J  fx)  =  A   cos  "$      (Bessel  function  of  1st  kind) 

Let  Y  (x)   =  A   sin  "$      (Bessel  function  of  2nd  kind) 
V  ^        V      —V 

Let  B  (t)  =  y7rt72  A  (t),  where  t  =  l/x  . 
Then  it  can  be  shown  that  near  t  =  0,  B  (t)  is 
approxiFiately  given  by  the  solution  of 

(1)  l^ii(il  +1   [i„(^2    1)^2]  B  (t)  -  [B  ft)]-^  =  0 

dt^      t^  ^       V         V 

or   (IM   t^  B^(t)  B;;(t)  +  [l-(v2-  l)t2]  Bj;(t)    -  t^  =  0 

and  that  there  is  a  solution  which  satisfies: 

■^v^^  2      4  l4 

(2)  — r —  =  1  +  a„t   +  aut   +  ...  +  a^ht   +  ...  . 
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It  can  also  be  shown  that  "$   satisfies 

(3)  -^t-  =     "  [B^(t)]-' 

or   (3M      B^(t)  ~^^   =  -  1 
•  V      dt 

and  that  there  is  a  solution  which  satisfies 

(4)  t[|^  +  (v+ I)  |]   =   1  +  b^t^  +  bi^t^  +  ... 

+  D-|  jfTj  +  •  .  .   • 

The  problem  was  to  find  the  coefficients   ap,...,a^2i 

2   1 
and  bp^  •  -  «  j^iu    i^^  terms  of  a  =  (v  =  -jr)  .   B  (t)  was  set 

3  IS 

equal  to  t  +  a^t   +  ...  +  a^ k t  ^  and  then  this  expression 

was  used  to  build  up  the  expression  on  the  left  side  of  (1') 

In  the  resulting  expression,  called  E,  the  coefficient  of 

the  lowest  power  of  t  was  set  equal  to  zero  and  solved 

for  ap.   E  was  then  adjusted  by  erasing  the  constant  term 

and  dividing  by  the  smallest  non  zero  power  of  t  appearing 

on  E.   The  solution  for  ap  was  then  substituted  into  E  and 

the  process  repeated  for  aw  instead  of  ap,  etc.   At  various 

stages  in  building  up  the  expression  (1')^  a  truncating 

routine  was  used  to  erase  terms  in  the  expression  which 

had  a  power  of  t  larger  than  the  smallest  power  of  t  which 

was  multiplied  by  a-,2,. 

In  order  to  obtain  the  coefficients  bp,  .  . .  ,b-,  ^ , 

equation  (4)  was  solved  for  ^  ,  differentiated  once  and 

2 
multiplied  by  t  ,  giving  a  pencil  and  paper  e:xpression: 

2  ^v  2  14 

[t^  -^]      =      ~   1   +   bpt^  +  .  .  .  +  13  b^i^t-"^ 

=  24  - 


which  was  set  up  as  a  computer  expression.   Then,  in  the 
computer,  the  expression 

(A)  [t^^]  -B^Ct)  +  t^  =  0 

which  is  equivalent  to  equation  (3')^  was  built  up.   After 
truncating  according  to  the  above  rule  with  respect  to  b-,  h, 
the  resulting  coefficients  of  t  on  (A)  were  solved  for 
bp,...,b^2,  in  a  manner  similar  to  that  used  for  obtaining 

ap,  •  •  •  ,  ^n  li  • 

The  coefficients  a^,  ...,a^h  and  bp,  ...,b,|.  are  given 

here  for  reference.   ap,...,a-,Q  and  bp,  ...,b^Q  were  published 

4  ?  1 

in  [11].   There,  the  coefficient  of  a   (with  a  =  {v   -  -^)) 

in  bo  was  printed  as  1/1792;  it  should  have  been  10/1792.. 

ap  =   1/4  a 

a^  =   5/32  a^   -  J>/8  a 

ag  =   15/128  a^   -  57/52  a^   +   I5/8  a 

ag  =   195/2048  a^  -  6II/256  a^  +  I82I/I28  a^  -   5I5/I6  a 

a^Q  =   665/8192  a^  -  4199/1024  a^  +  298II/5I2  a^   -  2225/8  a^ 

+  2855/8  a 
a^2  =   0O7O81604  a^   -   605262529  a^  +   17^.088l5  a^ 

-  1898.7920  a^   +  8025.559^  a^  -  9745.5125  a 
a^l^  =   .06522861  a^  -  9.0880585  a^  +   428.05176  a^ 

-  8787.6775  a^  +  82500.917  a^  -   525165.7^  a^ 

+  580067.19  a 
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b 


2 


=   1/2 


a 


=   1/24  a^  -  1/4 


a 


b 
b 


10 


12 


^14 


1/80  o?  -   7/20  a^  +  3/4  a 

10/1792  a^  -  95/224  a^  +  807/224  a^  -  315/56  a 
7/2304  a^  -  35/72  a^  +  1975/192  o?   -  58  a^  +  315/^ 
-.001864346  a   -  .54o66o49  a^  +  22.8405526  et^ 

-  303.317471  o?   +  1405.4829  a^  -  1771.8750  a 
.0012394837  oJ  -    .58999398  o^   +  43.573618  a^ 

-  1112.5907  a^  +  11695.116  a^  -  48514.975  a^ 
+  58471.875  a 


a 


In  the  expressions  for  a-,p,  ^iii^  ^ip^  ^i  ij.  -^^  -^^ 
probable  that  the  eighth  significant  figure  is  wrong 
because  single  precision  floating  point  arithnratic  was 
used. 


D.   A  Partial  Differential  Equation, 


(L(7r),l) 


u  =  0 


u  =  1  =  u(L(y),y) 


(L(0),0) 


Figure  1 


This  problem,  suggested  by  Professor  Herbert  Keller, 


is  to  solve  Laplace's  equation: 

(1.0)  1^  +^   =   0 


hyi 


^1 
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in  the  shaded  region  of  Figure  1.   The  solution  u(x,y)  is 
.subgectto  the  boundary  conditions: 


(1.1) 

u(0,y)      = 

0 

(1.2) 

u(L(y),y) 

= 

(1.3) 

Uy(x,0)      = 

0 

(1.4) 

u    (x,7r)      = 

0 

The  boundary  on  the  right  is  given  by  a  power 
series  in  £  with  Fourier  series  coefficients: 


(2.0)  X  =  L(y)   =  2   L.(y)  s 

j=0  ^ 


J 


--     (2,1)  •'-'o^'^^  ^      ^  (const.) 

00 

(2<,2)  L.(y)   =  Z  a.^   cos  ky 

J       k=0   J^ 

The  solution  of  eq.  (1.0)  will  be  obtained  by 
separation  of  variables,  then  the  problem  will  reduce 
to  determining  coefficients  of  series —  the  job  done 
by  the  algebraic  manipulator.   The  formulation  of  this 
solution  is  also  due  to  Professor  Keller. 

Let  u  =  X(x)*Y(y);   then  substituting  for  u  in 
(1.0)  leads  to 

ir         ti 

X  /X  +  Y  /Y  =   0 
and  then  since  the  first  term  is  a  function  of  x 
alone  and  the  second  term  a  function  of  y  alone; 

tt  p       "         p 

y  /Y  =   -  A       X  /X  =   A 

By  solving  these  ordinary  differential  equations,  the 
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solutions  are  found  to  be 

(3.1)  Y(y)   =  a  cos  Ay  +  p  sin  Ay 

(3.2)  X(x)  =      ^   Ae^^  +  Be"^^       A  ^  0 

A    +  Bx         A  =  0 

Because  of  (1.3),  It  Is  necessary  that  Y'(0)  =  0, 
which  from  (3'1)  Implies  that  p  =  0.   From  (1.4), 
Y'(7r)  =  0,  which  Implies  either  a  =  0  or  sin  Att  =  0 . 
If  a  =  0,  then  Y(y)  s  0  and  then  u  s  o  and  (1.2) 
couldn't  be  satisfied  by  u(x,y)  »   If  sin  Att  =  0,  then 
the  possible  values  of  A  are  0,  ±1,  ±2,     ...    .      Now,  we 
assume  that  the  solution  u(x,y)  to  (1.0)  with  the 
prescribed  boundary  conditions  Is  the  sum  of  all  functions 
u  -  X(x)-Y(y)   I.e., 

n=+^ 

u(x,y)   =   (Aq+BqX)  +  XI  (A^e""^  +  B^e""^)  cos  ny  . 

n=— CO 

n^^O 

Because  of  the  symmetric  way  in  which  the  Integer  n 

enters  the  expression  for  u(x,y),  it  Is  possible  to  relabel 

the  A  's  and  B  's  and  get: 
n        n 

00 

(3.3)  u(x,y)  =  (A  +B  x)  +  2Z  (A  e^+B  e"^-^)cos  ny 

n=l 

From  (1.1)  and  (3.3)  we  find 

00 

(3.4)  u(0,y)   =  0  =  Aq  +2Z  ^^n+^n'  °°®  ^^  * 

n=l 
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¥e  assume  that  this  series  converges  uniformly  on 

the  interval  [0,ir],    then  Aq  =  0  and  A^+   B  =  0  for  n  7^  0 , 

Then  (3.3)  becomes 

00 

(3.5)  u(x,y)   =  '^n^  +yZ   B^(e~^-e'^)  cos  ny  . 

^    n=l  ^ 

By  relabeling  the  B. 's,  this  becomes 

00 

(3.6)  u(x,y)   =   bp^x  +  >   b   sinh  nx  cos  ny  , 

n=l   ^ 

and  (3.6)  satisfies  (1.0),  (1.1),  (1.3)  and  (1.4). 

The  problem  is  now  reduced  to  finding  the  b  's  so 
that  (1.2)  is  satisfied  when  the  right  boundary  is  given 
by  the  equations  (2.0),  (2.1)  and  (2.2).   We  assume  that 


(3.7)  ^^     =     IZ  ^kn  ^         k=0,l,...  . 

^     j=0  '^^ 

The  problem  is  now  reduced  to  determining  the  constants 

b,  .  in  terms  of  the  constants  a  and  a  .1^.   Prom  (1,2),  (2.0), 

and  (3.6)  we  have  at  the  right  boundary 


(3.8)   u(L(y),y)  =  bQL(y)  +  ZZ  ^n  slnh(nL(y)  )  cos  ny-1 

n=  1 

or 


(3.9)   1   -  ^Q  ;^  Lj(y) 


j=0 


00  00 

J 


>   b   sinh[  (na+  n  >   L  .  (y)e"^ )  ]  cos  ny  . 
"^  ^  >1   ^ 


Now,  using  sinh(x+y)  =  sinh  x  cosh  y  +  cosh  x  sinh  y 
this  becomes 

-   29  - 


(4.0)      b^  y~  L.(y)    e-^   +7""  b    [slnh  na   cosh(nT~L  .  (y)e'^ ) 
^    J=0      -^  n=l     ^  j^   '^ 

00 

+  cosh  na-slnh(n  ^        L  .  (y)e'^ )  ]  cos   ny  -   1     =      0    , 

3=1      ^ 


where  the  quantity  n  y        L.(y)e'^   can  be  made  arbitrarily 
small  by  choosing  e  small,  in  which  case  the  power  series 
expansions  for  cosh  z  and  slnh  z  give  good  approximations. 

As  an  example  of  how  the  bi^-'s  are  obtained,  the 
equation  resulting  from  the  constant  term"  with  respect 
to  e  follows.   Prom  (4.0),  we  obtain 


(5.0)     0  =   -1  +  bp|L„(y)  +  y        b   slnh  na  cos  ny  , 

n=l  ^ 

which  becomes,  after  removing  terms  in  e  from  bp^,  b   and 


0'  "n 


Lo(y) 


(5.1)  0  =   -  1  +  bpiQa  +  y        slnh(na)  b  q'cos  ny  . 

n=l 

This  equation  is  to  hold  for  all  y  In  [0,Tr]   and  thus 
assuming  uniform  convergence  for  the  series  on  the  right: 

(5.2)  b^Q  =  1/a     and     b^^  =  0   for  n  =  1,2,, 

In  the  machine  computation,  only  two  terms  were  taken 
in  either  the  fourler  series  or  power  series  of  a  function. 
Specifically  the  following  substitutions  were  made  in  (4.0): 

^k  =   ^kO  "^  ^kl^  "^  ^k2^^     ^°^   ^  ^   0,1, 2j  bj^=0,  k  >  2, 
L  (y)   =   a  for  j  =  0 

J 

L.(y)    =  ct  .Q^  .-lCOs   y  +  aj2COS    2y   for    j   =    1, 2;Lj  (y)=0,  j>2. 
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CXD  .  2  ,    „ 

cosh  n  (  YZ   L,(y)  e^)      =      1  +  (n  >_  L  (y)e'')V2  , 

for  n  =  1,2;  and  =0  for  n  >2, 

00  .  2  . 

sinh  n  (  7^  L  (y)e'^)   =   (n  ^2  ^^(y)    e"^)  > 
3=1      ^  j=l   ^ 

for  n  =  1,2; and  =  0  for  n>2. 

After  the  new  expression  for  (4.0)  was  expanded  and 
terms  collected,  terms  containing  powers  of  e  greater  than 
e   were  dropped.   Terms  which  contained  products  of  cosines 
were  reduced  by  using  the  formula 

-  ■   —   2  cos  A  cos  B  =   cos  (A+B)  +  cos  (A-B) 

If  A+B  >  2y,  then  the  term  in  A+B  was  dropped.  After  all 
this,  we  were  left  with  an  expression  which  resulted  from 
(4.0)  and  which  was  of  the  form: 

Aq  +  A^e  +  Age^  =   0 
where 

Aq   was  of  the  form  Cqq  +  C^-j^cos  y  +  C^g  cos  2y 

A.   was  of  the  form  C^^   +  C-j^^cos  y  +  C^^  °°s  2y 

Ap  was  of  the  form  CpQ  +  C22_cos  y  +  C^^   ^^^    ^y 

where,  in  turn,  the  C.  .'s  were  linear  combinations  of  the 
b.  .'s.   These  relations  led  to  the  collection  of  equations 
to  be  solved  for  the  b^.'s: 
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(6.0)   Cqq  =  0,   Cq^  =  0,   Cq2  =  0.   C^Q  =  0,   c-L3_  =  0, 

Cnp  ~  ^>       '^20  "^       21  ~       22  "^    * 

The  unknowns  IDqq,  b^^,  Td^q^  ^01'  ''^11'  ^21^  ^02' 
b^p,  bpp  were  involved  linearly  in  the  expressions  for 

^00^  ^01'  ^02'  ^10'  ^11^  ^12-'  ^20'  ^21'  ^22°   These 
expressions  C  .  were  set  equal  to  zero  and  solved  for  b  . . 

in  the  order  given  above.   The  solution  for  each  b.  .  was 

substituted  into  each  of  the  expressions  following  C  . . , 

J 

and  the  process  continued. 

6 .   Conclusions . 

All  four  problems  discussed  in  Section  5  support  the 
conclusion  that  it  is  possible  to  use  computers  to  do 
algebraic  problems  which  would  otherwise  be  extremely 
tedious,  if  not  impossible.   The  first  problem  —  which 
is  a  larger  algebraic  problem  than  that  which  can  be 
reliably  done  by  hand  --  didn't  strain  the  capabilities  of 
this  primitive  algebraic  manipulator.   Thus  the  conclusion 
that  it  is  now  possible  to  do  algebraic  manipulation  of 
extraordinary  magnitude  is  justified. 

The  anticipated  difficulties  relating  to  factorization 
of  algebraic  expressions  never  appeared.   If  the  first 
problem  were  larger,  the  factored  coefficients  would  have 
been  introduced  into  the  machine  computation  for  expansion 
and  then  comparison.   Of  course,  the  ability  to  do  factoriza- 
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tion  would  be  useful. 

In  numerical  programming,  It  would  be  reasonable  to 
expect  an  experimental  set  of  routines  to  run  perhaps  twice 
as  long  as  a  polished  system.   However,  in  the  non-niimerical 
area,  especially  using  list  processors,  it  is  possible  to 
waste  factors  of  100  in  running  time.   The  basic  reason  for 
these  large  factors  in  inefficiency  is  that  the  number  of 
operations  in  non-numerical  algebraic  problems  tends  to  be 
described  by  combinatorial  factors  so  that  inefficiencies 
are  also  combinatorial  instead  of  linear.   The  use  of  a  list 
processor  essentially  changes  the  time  scale  by  a  factor  of 
about  10  -  50  because  operations  analogous  to  CLA  or  ST0 
on  the  7094  require  whole  subroutines  for  their  execution. 
These  two  facts  indicate  that  efficiencies  which  are 
introduced  into  such  systems  will  yield  dramatic  improve- 
ments . 

While  our  choice  of  SLIP  and  our  notation  were 
warranted  because  they  allowed  us  to  get  problems  written 
and  debugged  quickly,  we  see  that  for  practical  applica- 
tions in  the  future  we  will  need  to  require  the  routines 
using  a  more  compact  notation  and  a  faster  list  processor. 
Such  a  list  processor,  NU  SPEAK  [9],    developed  here,  seems 
eminently  suitable  for  this  task,  and  we  therefore  intend 
to  use  it  as  the  basis  for  a  new  version  of  SYMB0LANG . 
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